import pandas as pd
import numpy as np
import geopandas as gpd # Vector geospatial operations
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.warp
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")
/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow. warnings.warn(
['lds-nz-road-centrelines-topo-150k-FGDB.zip', 'subnational-population-projections-2018base-2048.xlsx', 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip', 'statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip', 'statsnzregional-council-2021-clipped-generalised-FGDB.zip', 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']
First, read regional council bounds. This geometry will be used to clip NZ-wide datasets to just the region of interest, Auckland
%%time
REGC = gpd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)
epsg:2193 CPU times: user 1.84 s, sys: 878 ms, total: 2.72 s Wall time: 14.1 s
Load the LRIS Land Cover Database (downloaded in GDB format from https://lris.scinfo.org.nz/layer/104400-lcdb-v50-land-cover-database-version-50-mainland-new-zealand/)
%%time
df = gpd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")
CPU times: user 1min 47s, sys: 2.71 s, total: 1min 50s Wall time: 1min 50s
print(df.columns)
print(df.crs)
display(df.sample(5))
Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
dtype='object')
epsg:2193
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | ... | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 313790 | Landslide | Landslide | Landslide | Landslide | Landslide | 12 | 12 | 12 | 12 | 12 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb2000005253 | MULTIPOLYGON (((1160923.451 4898204.158, 11609... |
| 214669 | Gravel or Rock | Gravel or Rock | Gravel or Rock | Gravel or Rock | Gravel or Rock | 16 | 16 | 16 | 16 | 16 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000014095 | MULTIPOLYGON (((1887044.793 5620048.149, 18870... |
| 388703 | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | 69 | 69 | 69 | 69 | 69 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000182838 | MULTIPOLYGON (((1960574.060 5784716.815, 19605... |
| 107539 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | 71 | 71 | 71 | 71 | 71 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000202128 | MULTIPOLYGON (((1945776.018 5602220.374, 19457... |
| 75253 | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | 54 | 54 | 54 | 54 | 54 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2004-06-30T00:00:00 | lcdb2000135128 | MULTIPOLYGON (((1439088.933 5045590.762, 14391... |
5 rows × 24 columns
%%time
df = gpd.clip(df, AKL)
CPU times: user 46.4 s, sys: 28.3 ms, total: 46.5 s Wall time: 46.5 s
df.sample(5)
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | ... | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 257236 | High Producing Exotic Grassland | High Producing Exotic Grassland | Exotic Forest | Exotic Forest | Exotic Forest | 40 | 40 | 71 | 71 | 71 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000408113 | POLYGON ((1776033.445 5897655.722, 1776027.801... |
| 429399 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | 71 | 71 | 71 | 71 | 71 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000217879 | POLYGON ((1779737.510 5905755.113, 1779733.434... |
| 281720 | Built-up Area (settlement) | Built-up Area (settlement) | Built-up Area (settlement) | Built-up Area (settlement) | Urban Parkland/Open Space | 1 | 1 | 1 | 1 | 2 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30T00:00:00 | lcdb1000003636 | POLYGON ((1746719.161 5913844.184, 1746732.995... |
| 365670 | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | 52 | 52 | 52 | 52 | 52 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2015-06-30T00:00:00 | lcdb1000456637 | POLYGON ((1731793.016 5938194.674, 1731764.232... |
| 296981 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | 71 | 71 | 71 | 71 | 71 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000409942 | POLYGON ((1756300.691 5879431.457, 1756287.981... |
5 rows × 24 columns
df.Name_2018.value_counts()
Exotic Forest 3981 Indigenous Forest 3673 Manuka and/or Kanuka 2282 Broadleaved Indigenous Hardwoods 1788 Built-up Area (settlement) 1350 High Producing Exotic Grassland 1326 Mangrove 1151 Urban Parkland/Open Space 1099 Estuarine Open Water 441 Orchard, Vineyard or Other Perennial Crop 436 Short-rotation Cropland 362 Lake or Pond 326 Herbaceous Saline Vegetation 303 Low Producing Grassland 291 Gorse and/or Broom 287 Forest - Harvested 266 Sand or Gravel 252 Deciduous Hardwoods 201 Surface Mine or Dump 132 Mixed Exotic Shrubland 120 Herbaceous Freshwater Vegetation 118 Transport Infrastructure 107 River 15 Flaxland 9 Gravel or Rock 9 Fernland 2 Matagouri or Grey Scrub 1 Name: Name_2018, dtype: int64
These classes are far too detailed - simplify to just Urban, Vegetation, Water, Other
def simplify_classes(code):
if code in [1, 2, 5]:
return 1, "Urban"
elif code in [68,69,71]:
return 2, "Vegetation"
elif code in [0,20,21,22,45,46]:
return 3, "Water"
else:
return 4, "Other"
summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
print(year)
class_year = f"Class_{year}"
df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
summary.append(df[class_year + "_simplified_name"].value_counts())
1996 2001 2008 2012 2018
pd.DataFrame(summary).plot.area()
<AxesSubplot:>
%%capture
# %%capture suppresses output
if not os.path.isfile("land_use.gif"):
ims = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
ax.set_title(year)
ax.figure.tight_layout()
canvas = ax.figure.canvas
canvas.draw() # draw the canvas, cache the renderer
image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(canvas.get_width_height()[::-1] + (3,))
ims.append(image)
imageio.mimsave("land_use.gif", ims, fps=1)
with open('land_use.gif','rb') as file:
display(Image(file.read()))
cols = [f"Class_{year}_simplified_code" for year in years]
cols
['Class_1996_simplified_code', 'Class_2001_simplified_code', 'Class_2008_simplified_code', 'Class_2012_simplified_code', 'Class_2018_simplified_code']
%%time
geocube = make_geocube(
vector_data=df,
output_crs="epsg:2193",
measurements=cols,
resolution=(-100, 100),
fill=0, # NaNs, like offshore areas, will be 0
)
geocube
CPU times: user 18.3 s, sys: 9.77 ms, total: 18.3 s Wall time: 18.3 s
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
spatial_ref int64 0
Data variables:
Class_1996_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2001_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2008_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2012_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2018_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])geocube.Class_2018_simplified_code.plot()
<matplotlib.collections.QuadMesh at 0x7faeb6597b20>
for year in years:
print(year)
outfile = f"output/land_use_{year}.tif"
if not os.path.isfile(outfile):
geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize
1996 2001 2008 2012 2018
%%time
TALB = gpd.read_file("input/statsnzterritorial-authority-local-board-2021-clipped-generalised-FGDB.zip!territorial-authority-local-board-2021-clipped-generalised.gdb")
TALB = gpd.clip(TALB, AKL)
TALB.head()
CPU times: user 1.2 s, sys: 9.79 ms, total: 1.21 s Wall time: 1.21 s
| TALB2021_V1_00 | TALB2021_V1_00_NAME | TALB2021_V1_00_NAME_ASCII | LAND_AREA_SQ_KM | AREA_SQ_KM | Shape_Length | geometry | |
|---|---|---|---|---|---|---|---|
| 2 | 00300 | Kaipara District | Kaipara District | 3108.960347 | 3108.960347 | 934117.838186 | MULTIPOLYGON (((1718576.593 5980699.195, 17185... |
| 4 | 01200 | Hauraki District | Hauraki District | 1270.133527 | 1270.133527 | 311282.487609 | POLYGON ((1801726.041 5899678.947, 1803788.474... |
| 5 | 01300 | Waikato District | Waikato District | 4403.797301 | 4451.122298 | 680627.805117 | POLYGON ((1747419.052 5870470.171, 1747413.783... |
| 44 | 07601 | Rodney Local Board Area | Rodney Local Board Area | 2275.114945 | 2275.114945 | 920471.582397 | MULTIPOLYGON (((1717309.547 5972989.294, 17173... |
| 45 | 07602 | Hibiscus and Bays Local Board Area | Hibiscus and Bays Local Board Area | 110.097301 | 110.097301 | 162589.255306 | MULTIPOLYGON (((1751233.745 5948474.046, 17512... |
# From https://www.stats.govt.nz/information-releases/subnational-population-projections-2018base2048#map
pop = pd.concat(pd.read_excel(
"input/subnational-population-projections-2018base-2048.xlsx",
sheet_name=["Table 5", "Table 6"],
skiprows=5,
usecols="A,B,G",
names=["area", "year", "population"],
nrows=231
)).fillna(method='ffill').reset_index()
pop.replace("Maungakiekie-Tamaki local board area", "Maungakiekie-Tāmaki local board area", inplace=True)
pop.year = pop.year.astype(str)
pop
| level_0 | level_1 | area | year | population | |
|---|---|---|---|---|---|
| 0 | Table 5 | 0 | Far North district | 1996 | 54500 |
| 1 | Table 5 | 1 | Far North district | 2001 | 56400 |
| 2 | Table 5 | 2 | Far North district | 2006 | 57500 |
| 3 | Table 5 | 3 | Far North district | 2013 | 60600 |
| 4 | Table 5 | 4 | Far North district | 2018 | 67900 |
| ... | ... | ... | ... | ... | ... |
| 457 | Table 6 | 226 | Franklin local board area | 2028 | 96900 |
| 458 | Table 6 | 227 | Franklin local board area | 2033 | 109400 |
| 459 | Table 6 | 228 | Franklin local board area | 2038 | 122000 |
| 460 | Table 6 | 229 | Franklin local board area | 2043 | 134500 |
| 461 | Table 6 | 230 | Franklin local board area | 2048 | 146900 |
462 rows × 5 columns
pop[pop.area.str.contains("local")].pivot(index='year', columns='area', values='population').plot()
<AxesSubplot:xlabel='year'>
pop = pop.pivot(index='area', columns='year', values='population')
pop
| year | 1996 | 2001 | 2006 | 2013 | 2018 | 2023 | 2028 | 2033 | 2038 | 2043 | 2048 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| area | |||||||||||
| Albert-Eden local board area | 86100 | 90000 | 96200 | 100000 | 103700 | 106900 | 111200 | 116500 | 121600 | 126300 | 130600 |
| Auckland | 1115800 | 1218300 | 1373000 | 1493200 | 1654800 | 1778700 | 1891800 | 2001800 | 2107000 | 2207800 | 2302900 |
| Devonport-Takapuna local board area | 51600 | 52300 | 55700 | 58500 | 60500 | 62700 | 65700 | 68800 | 71800 | 74800 | 77600 |
| Far North district | 54500 | 56400 | 57500 | 60600 | 67900 | 72400 | 75000 | 77300 | 79100 | 80600 | 81700 |
| Franklin local board area | 48400 | 52900 | 60600 | 68300 | 77700 | 85700 | 96900 | 109400 | 122000 | 134500 | 146900 |
| Gisborne district | 47200 | 45500 | 45900 | 47000 | 49500 | 51900 | 53100 | 54000 | 54600 | 55000 | 55200 |
| Great Barrier local board area | 1230 | 1100 | 940 | 950 | 960 | 1030 | 1080 | 1120 | 1150 | 1180 | 1200 |
| Hamilton city | 113500 | 121200 | 134800 | 150200 | 168600 | 183000 | 194400 | 205400 | 216000 | 226500 | 236600 |
| Hauraki district | 18550 | 18000 | 18300 | 18600 | 20700 | 21800 | 22200 | 22400 | 22300 | 22100 | 21800 |
| Henderson-Massey local board area | 82200 | 90900 | 103700 | 113500 | 124600 | 134100 | 142500 | 150500 | 158200 | 165700 | 172900 |
| Hibiscus and Bays local board area | 68000 | 75100 | 85200 | 94000 | 108500 | 117700 | 123500 | 126500 | 128900 | 130900 | 132400 |
| Howick local board area | 80900 | 97300 | 119100 | 135000 | 149400 | 160500 | 167800 | 174400 | 180400 | 185900 | 190900 |
| Kaipara district | 17800 | 17950 | 18550 | 20500 | 23700 | 25900 | 27200 | 28300 | 29100 | 29800 | 30300 |
| Kaipātiki local board area\n | 74700 | 78200 | 83600 | 87000 | 92900 | 97300 | 99400 | 101100 | 102500 | 103600 | 104700 |
| Kawerau district | 8120 | 7290 | 7150 | 6650 | 7460 | 7910 | 8000 | 8020 | 7970 | 7860 | 7720 |
| Manurewa local board area | 60100 | 69100 | 81600 | 87000 | 100900 | 109600 | 114000 | 118000 | 121700 | 124900 | 127600 |
| Matamata-Piako district | 30300 | 30300 | 31200 | 32900 | 35300 | 37000 | 37900 | 38700 | 39200 | 39500 | 39600 |
| Maungakiekie-Tāmaki local board area | 61700 | 66400 | 70400 | 73700 | 80500 | 86600 | 93800 | 101300 | 108800 | 116100 | 123000 |
| Māngere-Ōtāhuhu local board area\n | 61400 | 64600 | 72500 | 75300 | 82700 | 88700 | 92500 | 96300 | 99800 | 103100 | 105900 |
| Papakura local board area | 38000 | 38900 | 43100 | 48200 | 61100 | 70000 | 77000 | 81700 | 86100 | 90500 | 95000 |
| Puketāpapa local board area\n | 46100 | 49800 | 53900 | 56300 | 60900 | 64600 | 68800 | 73000 | 77100 | 81000 | 84600 |
| Rodney local board area | 39000 | 44100 | 51000 | 57300 | 69100 | 77600 | 87600 | 99600 | 111700 | 123800 | 135800 |
| Rotorua district | 66600 | 66900 | 68100 | 68400 | 74800 | 78900 | 80700 | 82200 | 83400 | 84200 | 84800 |
| South Waikato district | 25800 | 24200 | 23200 | 23200 | 24800 | 25800 | 26400 | 26800 | 27000 | 27100 | 27100 |
| Taupō district\n | 31600 | 32500 | 33400 | 34800 | 38600 | 40900 | 42000 | 42800 | 43300 | 43600 | 43800 |
| Tauranga city | 79900 | 93600 | 107000 | 119900 | 142100 | 156900 | 166300 | 175000 | 183300 | 191400 | 199100 |
| Thames-Coromandel district | 25400 | 25800 | 26700 | 27300 | 30700 | 32400 | 33100 | 33500 | 33500 | 33300 | 32800 |
| Upper Harbour local board area | 24300 | 33700 | 45100 | 56800 | 66800 | 75900 | 86000 | 95900 | 105800 | 115400 | 124900 |
| Waiheke local board area | 6680 | 7560 | 8150 | 8630 | 9360 | 9830 | 10250 | 10650 | 10900 | 11100 | 11200 |
| Waikato district | 52000 | 53700 | 59500 | 66500 | 78200 | 86100 | 92800 | 99300 | 105700 | 111900 | 117700 |
| Waipa district | 38400 | 40000 | 43700 | 48700 | 55000 | 59300 | 62100 | 64600 | 66900 | 68900 | 70700 |
| Waitematā local board area\n | 44900 | 52500 | 66600 | 81300 | 88500 | 94700 | 102800 | 111900 | 120700 | 129200 | 137400 |
| Waitomo district | 10000 | 9780 | 9680 | 9340 | 9630 | 9740 | 9730 | 9670 | 9540 | 9330 | 9070 |
| Waitākere Ranges local board area\n | 42100 | 44200 | 47500 | 50700 | 54200 | 56900 | 58300 | 59500 | 60300 | 60900 | 61100 |
| Western Bay of Plenty district | 35600 | 38900 | 42900 | 45400 | 53300 | 58100 | 60900 | 63300 | 65200 | 66700 | 68000 |
| Whakatāne district\n | 34200 | 34100 | 34500 | 34200 | 37100 | 38800 | 39300 | 39500 | 39500 | 39300 | 38900 |
| Whangārei district\n | 68400 | 70000 | 76500 | 83700 | 94100 | 101000 | 105500 | 109500 | 113100 | 116400 | 119300 |
| Whau local board area | 59700 | 66200 | 73100 | 76700 | 84100 | 89500 | 95000 | 100500 | 105900 | 111000 | 115900 |
| Ōpōtiki district\n | 9620 | 9490 | 9200 | 8780 | 9670 | 10250 | 10350 | 10400 | 10300 | 10150 | 9910 |
| Ōrākei local board area\n | 70900 | 73000 | 78400 | 83700 | 87700 | 91300 | 96200 | 101100 | 105800 | 110300 | 114500 |
| Ōtara-Papatoetoe local board area\n | 67800 | 70400 | 76600 | 80300 | 90500 | 97700 | 101500 | 104000 | 106000 | 107500 | 108600 |
| Ōtorohanga district\n | 9960 | 9590 | 9310 | 9590 | 10500 | 11050 | 11300 | 11550 | 11750 | 11900 | 12000 |
cols = pop.columns.tolist()
cols
['1996', '2001', '2006', '2013', '2018', '2023', '2028', '2033', '2038', '2043', '2048']
pop = pd.merge(TALB, pop, left_on=TALB.TALB2021_V1_00_NAME.str.lower().str.strip(), right_on=pop.index.str.lower().str.strip(), how="left")
pop[["TALB2021_V1_00_NAME"] + cols]
| TALB2021_V1_00_NAME | 1996 | 2001 | 2006 | 2013 | 2018 | 2023 | 2028 | 2033 | 2038 | 2043 | 2048 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Kaipara District | 17800 | 17950 | 18550 | 20500 | 23700 | 25900 | 27200 | 28300 | 29100 | 29800 | 30300 |
| 1 | Hauraki District | 18550 | 18000 | 18300 | 18600 | 20700 | 21800 | 22200 | 22400 | 22300 | 22100 | 21800 |
| 2 | Waikato District | 52000 | 53700 | 59500 | 66500 | 78200 | 86100 | 92800 | 99300 | 105700 | 111900 | 117700 |
| 3 | Rodney Local Board Area | 39000 | 44100 | 51000 | 57300 | 69100 | 77600 | 87600 | 99600 | 111700 | 123800 | 135800 |
| 4 | Hibiscus and Bays Local Board Area | 68000 | 75100 | 85200 | 94000 | 108500 | 117700 | 123500 | 126500 | 128900 | 130900 | 132400 |
| 5 | Upper Harbour Local Board Area | 24300 | 33700 | 45100 | 56800 | 66800 | 75900 | 86000 | 95900 | 105800 | 115400 | 124900 |
| 6 | Kaipātiki Local Board Area | 74700 | 78200 | 83600 | 87000 | 92900 | 97300 | 99400 | 101100 | 102500 | 103600 | 104700 |
| 7 | Devonport-Takapuna Local Board Area | 51600 | 52300 | 55700 | 58500 | 60500 | 62700 | 65700 | 68800 | 71800 | 74800 | 77600 |
| 8 | Henderson-Massey Local Board Area | 82200 | 90900 | 103700 | 113500 | 124600 | 134100 | 142500 | 150500 | 158200 | 165700 | 172900 |
| 9 | Waitākere Ranges Local Board Area | 42100 | 44200 | 47500 | 50700 | 54200 | 56900 | 58300 | 59500 | 60300 | 60900 | 61100 |
| 10 | Waitematā Local Board Area | 44900 | 52500 | 66600 | 81300 | 88500 | 94700 | 102800 | 111900 | 120700 | 129200 | 137400 |
| 11 | Whau Local Board Area | 59700 | 66200 | 73100 | 76700 | 84100 | 89500 | 95000 | 100500 | 105900 | 111000 | 115900 |
| 12 | Albert-Eden Local Board Area | 86100 | 90000 | 96200 | 100000 | 103700 | 106900 | 111200 | 116500 | 121600 | 126300 | 130600 |
| 13 | Puketāpapa Local Board Area | 46100 | 49800 | 53900 | 56300 | 60900 | 64600 | 68800 | 73000 | 77100 | 81000 | 84600 |
| 14 | Ōrākei Local Board Area | 70900 | 73000 | 78400 | 83700 | 87700 | 91300 | 96200 | 101100 | 105800 | 110300 | 114500 |
| 15 | Maungakiekie-Tāmaki Local Board Area | 61700 | 66400 | 70400 | 73700 | 80500 | 86600 | 93800 | 101300 | 108800 | 116100 | 123000 |
| 16 | Howick Local Board Area | 80900 | 97300 | 119100 | 135000 | 149400 | 160500 | 167800 | 174400 | 180400 | 185900 | 190900 |
| 17 | Māngere-Ōtāhuhu Local Board Area | 61400 | 64600 | 72500 | 75300 | 82700 | 88700 | 92500 | 96300 | 99800 | 103100 | 105900 |
| 18 | Ōtara-Papatoetoe Local Board Area | 67800 | 70400 | 76600 | 80300 | 90500 | 97700 | 101500 | 104000 | 106000 | 107500 | 108600 |
| 19 | Manurewa Local Board Area | 60100 | 69100 | 81600 | 87000 | 100900 | 109600 | 114000 | 118000 | 121700 | 124900 | 127600 |
| 20 | Papakura Local Board Area | 38000 | 38900 | 43100 | 48200 | 61100 | 70000 | 77000 | 81700 | 86100 | 90500 | 95000 |
| 21 | Franklin Local Board Area | 48400 | 52900 | 60600 | 68300 | 77700 | 85700 | 96900 | 109400 | 122000 | 134500 | 146900 |
d = pop[cols].to_numpy()
d.min(), d.mean(), d.max()
(17800, 84495.2479338843, 190900)
%%capture
# %%capture suppresses output
if not os.path.isfile("pop.gif"):
ims = []
for col in cols:
ax = pop.plot(column=col, legend=True, vmin=0, vmax=d.max(), cmap='OrRd')
ax.set_title(col)
ax.figure.tight_layout()
canvas = ax.figure.canvas
canvas.draw() # draw the canvas, cache the renderer
image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(canvas.get_width_height()[::-1] + (3,))
ims.append(image)
imageio.mimsave("pop.gif", ims, fps=1)
with open('pop.gif','rb') as file:
display(Image(file.read()))
%%time
popcube = make_geocube(
vector_data=pop,
measurements=cols,
like=geocube, # Ensures dimensions match
fill=0 # NaNs, like offshore areas, will be 0
)
popcube
CPU times: user 1.88 s, sys: 0 ns, total: 1.88 s Wall time: 1.88 s
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
1996 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2001 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2006 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2013 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2018 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2023 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2028 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2033 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2038 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2043 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
2048 (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])for col in cols:
outfile = f"output/pop_{col}.tif"
if not os.path.isfile(outfile):
# byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
popcube[col].rio.to_raster(outfile, dtype=np.uint32)
%%time
roads = gpd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")
CPU times: user 11.8 s, sys: 10.8 ms, total: 11.8 s Wall time: 11.8 s
%%time
akl_roads = gpd.clip(roads, AKL)
CPU times: user 23.1 s, sys: 0 ns, total: 23.1 s Wall time: 23.1 s
# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway
| t50_fid | name_ascii | macronated | name | hway_num | rna_sufi | lane_count | way_count | status | surface | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 512 | 100120610 | KAIPARA COAST HIGHWAY | N | KAIPARA COAST HIGHWAY | 16 | 3007739 | 2 | None | None | sealed | LINESTRING (1732000.000 5944172.070, 1732048.5... |
| 2933 | 3198057 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 2 | None | None | sealed | LINESTRING (1748581.508 5968975.145, 1748558.4... |
| 2934 | 3198059 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 2 | None | None | sealed | LINESTRING (1748171.047 5971284.152, 1748129.9... |
| 3320 | 3200754 | PAERATA ROAD | N | PAERATA ROAD | 22 | 3000260 | 2 | None | None | sealed | LINESTRING (1767236.112 5888088.508, 1767244.3... |
| 3324 | 3200792 | UPPER HARBOUR MOTORWAY | N | UPPER HARBOUR MOTORWAY | 18 | 3047073 | 4 | None | None | sealed | LINESTRING (1747954.314 5927269.837, 1747970.0... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 138240 | 100048291 | AUCKLAND-WAIWERA MOTORWAY | N | AUCKLAND-WAIWERA MOTORWAY | 1 | 3067966 | 7 | None | None | sealed | LINESTRING (1755881.018 5922863.734, 1755886.4... |
| 138301 | 100048432 | AUCKLAND-HAMILTON MOTORWAY | N | AUCKLAND-HAMILTON MOTORWAY | 1 | 3017109 | 1 | None | None | sealed | LINESTRING (1765115.647 5909916.697, 1765092.7... |
| 138337 | 100048532 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 4 | None | None | sealed | LINESTRING (1748892.089 5949596.727, 1748892.0... |
| 138369 | 100048589 | PORT ALBERT ROAD | N | PORT ALBERT ROAD | 16 | 3013274 | 2 | None | None | sealed | LINESTRING (1734173.019 5980575.187, 1734175.6... |
| 138680 | 100118365 | SOUTH-WESTERN MOTORWAY | N | SOUTH-WESTERN MOTORWAY | 20 | 3018532 | 4 | None | None | sealed | LINESTRING (1760066.252 5908184.133, 1760043.9... |
426 rows × 11 columns
mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()
mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()
ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)
%%time
mway_cube = make_geocube(
vector_data=mway,
measurements=["lane_count"],
like=geocube, # Ensures dimensions match
fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube
CPU times: user 203 ms, sys: 162 µs, total: 203 ms Wall time: 200 ms
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
lane_count (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)
src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
src_ds.RasterXSize, src_ds.RasterYSize, 1,
gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")
(1320, 1001)
Text(0.5, 1.0, 'Distance (m)')
%%time
skytower = [-36.84838748948485, 174.7621736911587]
skytower = gpd.points_from_xy(x=[skytower[1]], y=[skytower[0]])
cbd = gpd.GeoDataFrame([{"name": "Skytower", "value": 1}], geometry=skytower, crs="EPSG:4326").to_crs(epsg=2193)
display(cbd)
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
vector_data=cbd,
like=geocube, # Ensures dimensions match
fill=0
)
display(cbd_cube)
outfile = "output/cbd.tif"
if not os.path.isfile(outfile):
cbd_cube.value.rio.to_raster(outfile, dtype=np.byte)
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
src_ds.RasterXSize, src_ds.RasterYSize, 1,
gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")
| name | value | geometry | |
|---|---|---|---|
| 0 | Skytower | 1 | POINT (1757109.057 5920497.145) |
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
value (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])(1320, 1001) CPU times: user 360 ms, sys: 10 ms, total: 370 ms Wall time: 365 ms
Text(0.5, 1.0, 'Distance (m)')
bounds = AKL.total_bounds.tolist()
bounds
[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]
zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles
['EJ.tif', 'DM.tif', 'EL.tif', 'DL.tif', 'DJ.tif', 'FK.tif', 'DK.tif', 'EK.tif', 'FL.tif']
tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transform = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)
print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')
-32767 -18316.517014943143 697 (1, 1320, 1018)
<AxesSubplot:>
cbd_dist = rasterio.open("output/cbd_dist.tif")
cbd_dist.read().shape
(1, 1320, 1001)
NODATA = -32767
warped_DEM = np.full_like(cbd_dist.read(), NODATA, dtype=np.int16)
rasterio.warp.reproject(DEM,
warped_DEM,
src_transform=transform,
src_nodata=NODATA,
dst_nodata=NODATA,
src_crs=tile_datasets[0].crs,
dst_crs=cbd_dist.crs,
dst_transform=cbd_dist.transform
)
print(DEM.shape, warped_DEM.shape)
(1, 1320, 1018) (1, 1320, 1001)
meta = tile_datasets[0].meta
print(meta)
meta.update({
"dtype": "int16",
"height": cbd_dist.height,
"width": cbd_dist.width,
"transform": cbd_dist.transform
})
print(meta)
outfile = "output/dem.tif"
if not os.path.isfile(outfile):
with rasterio.open(outfile, "w", **meta) as dest:
dest.write(warped_DEM)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703900.0,
0.0, -100.0, 6002400.0)}
gdal.DEMProcessing('output/slope.tif', outfile, 'slope')
slope = rasterio.open("output/slope.tif").read()
print(np.nanmin(slope), np.nanmean(slope), np.nanmax(slope), slope.shape)
rasterio.plot.show(np.where(slope>=0, slope, np.nan), cmap='terrain')
-9999.0 -5778.3354 52.492313 (1, 1320, 1001)
<AxesSubplot:>
rasters = [rasterio.open(f"output/{f}") for f in os.listdir("output") if f.endswith(".tif")]
params = set([(r.shape, r.res, r.crs, r.count, r.bounds) for r in rasters])
print(params)
# Assert all rasters have the same shape, pixel size, CRS, number of bands, and bounds
assert len(params) == 1
{((1320, 1001), (100.0, 100.0), CRS.from_epsg(2193), 1, BoundingBox(left=1703900.0, bottom=5870400.0, right=1804000.0, top=6002400.0))}
!ls -Ggh output
total 66M -rw-r--r-- 1 1.3M Apr 19 14:35 cbd.tif -rw-r--r-- 1 2.6M Apr 19 14:35 cbd_dist.tif -rw-r--r-- 1 2.6M Apr 19 11:42 dem.tif -rw-r--r-- 1 1.3M Apr 15 14:25 land_use_1996.tif -rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2001.tif -rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2008.tif -rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2012.tif -rw-r--r-- 1 1.3M Apr 15 14:25 land_use_2018.tif -rw-r--r-- 1 1.3M Apr 15 14:37 mway.tif -rw-r--r-- 1 2.6M Apr 19 14:18 mway_dist.tif -rw-r--r-- 1 5.1M Apr 19 14:05 pop_1996.tif -rw-r--r-- 1 5.1M Apr 19 14:10 pop_2001.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2006.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2013.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2018.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2023.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2028.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2033.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2038.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2043.tif -rw-r--r-- 1 5.1M Apr 19 14:11 pop_2048.tif -rw-r--r-- 1 5.1M Apr 19 14:18 slope.tif